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Abstract 

Recently a path integral formalism has been proposed by the anthor which gives the time evo¬ 
lution of moments of slow variables in a Hamiltonian statistical system. This closure relies on 
evaluating the informational discrepancy of a time sequence (path) of approximating densities from 
the Liouvillian evolution that an exact density must follow. The discrepancy is then used to weight 
all possible paths using a generalized Boltzmann principle. That formalism is extended here to 
deal with more general and realistic autonomous dynamical systems. There the divergence of the 
time derivative of dynamical variables need not vanish as it does in the Hamiltonian case and this 
property complicates the closure derivation. Many interesting and realistic applications are covered 
by this new formalism including those describing realistic turbulence and the relevant specifics of 
this situation are outlined. The practical issues associated with the implementation of the outlined 
formalism are also discussed. 
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I. INTRODUCTION 


The statistical closure of the dynamical systems underlying turbulent phenomenon is a 
long standing and practically important problem in mathematical physics. Recently the 
author has proposed a new formalism to deal with Hamiltonian systems using a Wiener 
path integral 1|. This formalism builds on earlier work by Turkington Q which has been 
validated numerically for first moments using an inviscid truncated Burgers system j^. The 
long term aim of this work is to apply the closure methodology to realistic turbulent system 
of various kinds. Consequently it would be desirable if the Hamiltonian restriction could be 
relaxed in order to deal with systems with dissipation and external forcing. 

As is well known (see e.g. Chapter 7) a large number of important conservative fluid 
dynamical systems can be cast into Hamiltonian form. Depending on the manner in which 
this is done and the particular fluid system, these may be either canonical or non-canonical. 
Specifically this means that as dynamical systems, they may be written as Q 


{Q,P] = J. 


dP dQ 


dxj dxi 


where H is the Hamiltonian and the Poisson bracket {,} is antisymmetric and satisfies 
the (Lie) Jacobi identity. Due to the bracket properties, the covariant tensor Jij is anti- 


:ic j^. 


symmetric [6|J. Furthermore an appropriate non-singular transformation of dynamical vari¬ 
ables allows it to be brought into the form of a constant anti-symmetric matrix (see ^ 
p333-337). For such a set of variables it is trivial to confirm that 

dA 


dxi 


= 0 


( 1 ) 


and since Ai is a covariant vector, this non-divergent condition holds for any non-degenerate 
choice of variables. As is also very well known such a condition is crucial for establishing that 
the statistical density evolution equation for an ensemble of dynamical system trajectories 
is the Liouville equation; 

Qt + Ai-^g = 0 ( 2 ) 

This equation implies easily that any function of the density also satisfies the same equa¬ 
tion. Furthermore condition ([1]) implies that the Liouvillian operator L = is anti- 

Hermitean with respect to the real inner product (see below for more detail on the more 
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Q 


general case). Both of these properties were needed to derive the closure in []| and so the 
non-divergence condition ([1]) was obviously essential. 

Consider now a conventional turbulent fluid system. As noted in ^ p233, this may often 
be described by the following dynamical system 


dxi 

dt 


Ai{x) + a{i)xi + F{i) 


(3) 


where brackets indicate no summation convention and where the A* belong to a conservative 
system which usually can be brought into Hamiltonian form as discussed above. The second 
term on the right denotes a (linear) dissipation (i.e. a < 0) while the third denotes a possible 
external and prescribed injection of energy. Such a system is often called a forced dissipative 
model of turbulence. If F is absent we obtain a decaying turbulence model. The differences 
of these two situations in the context of the present closure are discussed further below. The 
Navier Stokes equation, appropriately truncated spectrally, may be cast into the form of ([2]) 
with the coefficients a{i) increasing strongly with increasing wavenumber i. Physically the 
dissipation there is intended to represent molecular diffusion. Other turbulent systems such 
as a quasi-geostrophic fluid may have quite different forms for both the dissipation and Aj. 
In anycase the generic form of the dynamical system equations implies that 


dAj 

dxi 


< 0 


which must be true no matter what dynamical variables are chosen to describe the system. 
As discussed in detail below this non-zero divergence means that a non-trivial function 
of the statistical density no longer satisfies a generalized Liouville equation and also that 
the generalized Liouvillian operator is no longer anti-Hermitean. We therefore extend our 
closure derivation to an autonomous rather than a Hamiltonian dynamical system to deal 
with this more realistic situation. Our closure will be seen to be appropriate for a completely 
general autonomous dynamical system not just that of (jS]) and so further applications beyond 
conventional turbulence systems will also be possible. 

In section 2 we briefly review the earlier work mainly from a motivational and scene 
setting perspective before providing a complete derivation of our formalism for a general 
autonomous dynamical system in Sections 3 and 4. A discussion of applications may be 
found in Section 5. The reader is referred to the earlier papers mentioned above for more 
details on the Hamiltonian system closure methodology. 
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II. CONCEPTUAL DESCRIPTION OF THE CLOSURE FORMALISM 


A. Onsager-Machlup Theory 


It has been well known since early in the history of statistical physics that the probability 
of a fluctuation from an equibrium state at a fixed time is given by the Boltzmann principle 

p(A) = (7exp(-(7(A)) 

where a is the entropy function for the system. Onsager and Machlup (OM) Q general¬ 
ized this principle in the early 1950s and dehned a positive likelihood weight W on temporal 
paths X{t) of near equilibrium states 

r*t 


W [A(t)] = exp J ^ [A(f)] dt^ 


(4) 


The “Lagrangian” was assumed to be of the form; 

^[Xit)] = l[x-UXyg[x-UX 

where U and g are constant matrices with the latter importantly assumed to be positive 
dehnite. 

It is immediately clear from the latter assumption that the path of greatest weight is 
simply given by 

X = UX (5) 

which is, in general, a relaxation to equilibrium for appropriate choices for U. Furthermore 
the probability of a particular A at time T is obtained by integrating over the likelihood 
weights of all paths ending at A(T): 


p{X (T)) = C / DXW [A] 


( 6 ) 


which is a functional or path integral over all possible paths which end at A (T). As 
is well known 8j this situation corresponds with an Ornstein Uhlenbeck stochastic process 
and the most likely A also follows equation (|5]). Such a trajectory is called thermodynamical 
since it is most probable. 

Many attempts have been made to formulate far from equilibrium versions of OM theory. 
The approach here follows three principles: 
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1. The Lagrangian ^ should be derivable from a fundamental information theoretic ar¬ 
gument relating to paths since it effectively generalizes the entropy function underlying 
the original Boltzmann principle. 

2. The precise form of the Lagrangian should be deducible from first principles and reflect 
the statistical properties of the slow or coarse grained variables for the statistical 
system. 


3. In appropriate limits thermodynamical relations of the type proposed recently by 


Ottinger 


^ should be recovered since such theories work well in practical applications. 


B. Trial density manifold 


Zubarev proposed the concept of a trial density as a way of approximating densities of 
non-equilibrium states. One specifies a set of moments of slow variables as thermodynamical 
variables and then uses a maximum entropy principle to associate these with a “trial” density. 
In general such trial densities q are not exact densities q for the given non-equilibrium 
state with the given moments. This situation contrasts with the equilibrium Gibbs density. 
The philisophical approach taken here therefore differs from that adopted in equilibrium 
statistical physics in that it is tacitly assumed that the exact density for a system will never 
generally be available and thus one must deal always with an approximation manifold of 
densities. On such a manifold the “veracity” of each density will be assigned a non-negative 
weight which we shall term a consistency distribution. It is analogous to a quantum 
wavefunction as we shall see below. 

The set of all trial densities with a prescribed set of moments can be given a Riemannian 
manifold structure by using the associated Fisher information matrix as a metric tensor. 
This manifold structure is the essence of the subject of information geometry ll| and shall 
become important to our discussion later. 
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III. APPLICATION TO AN AUTONOMOUS DYNAMICAL SYSTEM 


A. Liouvillian discrepancy 


The equation for the density evolution within a general autonomous dynamical equation 
is simply the Fokker-Planck equation without a diffusion term 12| i.e. 


Qt + Lg — 0 


OXi 


with summation convention and where the original dynamical system is 


dxi 

dt 


Ai{x) 


(7) 


As usual the formal adjoint of the operator L is given by 


L* = 



From this we deduce that 

RA- 

L + = (8) 

UXi 

Consider now a temporal trajectory or path through the trial density manifold. In general 
this will not satisfy equation ([7|). A measure of the discrepancy of a given path from such an 
evolution can be dehned as follows: Suppose the trial path p(f) is evolved forward in time 
an additional At according to the above equation resulting in 


Q'{t + At) = exp (—AtL) g{t) 



exp (AtT) g{t) 


According to elementary information theory the “distance” between this evolved density 
and the assumed trial density g{t + At) can be measured by their relative entropy i.e. by 
D {g\t + Af)||p(f + At)). 

Now we require the time evolution of I = log g to calculate the information loss rate 
dehned using the relative entropy. Unlike the Hamiltionian case for which = q arbitrary 
functions of the density do not satisfy equation 0 so we need to compute the evolution of 
I explictly using a Taylor series. 
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A density evolving according to ([7]) will satisfy 

(T -L*)g = -L^Q 

Farther from the form of L* it follows easily nsing the chain rnle that 

,dAi 


(T - L*) F{g) = F'{T- L*) p = -F'^g 


and so 


or 


(T -L*)l = 


Tl = L*l - 


dxi 

dAj 

dxi 


(9) 


now assnming that A does not depend on time explicitly (the antonomons assnmption) 
then we have after iteration 

p)A- 

T^l = {L*Y I - * 


dxi 


and so 




where we are expanding only to second order since this will be appropriate below. 


B. Evolution of a Liouville residual 


Consider now a trial density g which does not evolve according to eqnation ([7]) since it is 
constrained to he within the trial density manifold. Denote by angle brackets the expectation 
with respect to this density at a particnlar time. For a general random variable F we have 

d{F) 


dt 


-{LF) = {Ft) + j F{T-L*)g 

= {Ft) + j F{T-L*)i 
= (F,) + {FR') 

R' = {T-L*)i 


( 10 ) 


where we are nsing integration by parts to re-express L as the adjoint operator on the 
first line and nsing the chain rnle on the second line. Choosing F = 1 we obtain 




( 11 ) 
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where Ris a generalized Liouville residual since it vanishes if a density evolves according 
to o (see equation ([9])). We can rewrite (ITOl) as 

- (LF) = {Ft) + {FR) - (f^^' 


dt 

and now set F = R obtaining 


dxi 


- {LR) = {TR) + {R^) -(r 


so 


-{{L + T)R) = (R^) -{R 


dxi 

dxi 


or using 


-{{T-L*)R) = (R^) 
or using the definition of R and R' as 


(T - L*)R) = (i?2) + ( (T - L*) 


dxi 


( 12 ) 


C. Information loss 


Q) 


For one timestep At this can be expressed as (see 

IL = D{g{t + At)||^(t + At)) 

= J g{t + At) (^l{t + At) — l{t + At) 

where g and I are the (generalized) Liouville evolved density and log density over the 
time step with the starting density and log density being those drawn from the trial density 
manifold at that earlier time i.e. g{t) = g{t) l{t) = l{t). 

Using ([7]) and the log density counterpart ([9]) this may be re-expressed to second order 
as 

'dAi 


IL = 


^-AtL 


m I e-^^^*i{t) - I At^ + 


+ ... -e 


AtT 


lit) 


dxi 2 V dxi 

Using integration by parts the first operator in the integrand can be shifted to the right 
and converted to an exponential of a multiple of the adjoint operator so we get 


IL = 


-AtL* 


- I At^ + ^L* 


= ( Kt) - e 


-AtL* 


dxi 

, dAt (At? , 


dxi 


dxi 

+ ... 


+ ...] -e 


AtT 


m 


^At{T-L* 


'((«) 
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where the angle bracket is an expectation with respect to the trial density at time t and 
on the second line we are using the fact that T and L* commute due to the autonomous 
dynamical system assumption. Expanding out the remaining exponential operators as Taylor 


series and retaining only terms to second order: 


IL = ( l{t) - At^ - 

' ^ ^ dxi 2 

-l{t) - At{T - L*)i{t) - 


dAi \ . ^ N 2 ^ / dAi 

' + {Aty L*' 


dx 




dXi 


-{T-L*yi{t) 


= ( -mA + ALl. 


{Aty 


dxi 


{W) + oyAtf) 


dAi\ , dAi (Aty 

' + At^ + ^^ 


dxi 


dxi 


R^ + {T- L*) 


dA 

dxi 


where on the second line we are using equations fllip and (IT^ as well as the autonomous 
assumption. There is a remarkable number of cancellations and a very simple result which 


nicely extends the Hamiltonian case. Recall that R vanishes for time evolution according to 
the generalized Liouville equation ([7j). 


IV. PATH INTEGRAL FORMALISM 


A. The Lagrangian 


Since IL has a straightforward information theoretic interpretation we set therefore 




e(t) 


(13) 


and by analogy with the Onsager and Machlup approach, identify the action S for this 
path 

S=yJ 

as the analog of the entropy for our proposed path Boltzmann principle. To make further 
progress we need to specify more concretely the trial density. A convenient choice here is to 
use maximum entropy theory with respect to an equilibrium density obtaining 


p (A, x) = Z ^ (A) exp (A*(5 (x) — /3y{x)) 


(14) 


where the function ip is the analog of invariant quantities for a Hamiltonian system (often 
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energy in that case). It may only be known approximately as a quadratic function and need 
not be assumed neccessarily to be an invariant of the dynamical system. 

The vector Q is an appropriately chosen set of slow variables for the system whose mo¬ 
ments interest us. For practical reasons they are usually chosen to be quadratic functions 
since calculation of expectations with respect to a trial density often requires a Gaussian 
density for tractability. In the case of a forced dissipative turbulent system a non-trivial 
equilibrium density occurs and this may be used empirically (via a direct numerical simula¬ 
tion) to derive '0(x). In the case of decaying turbulence when the forcing is absent then one 
expects the asymptotic state to be one of no motion meaning that the Lagrange multipliers 
A associated with quadratic variables might be expected to diverge as time proceeds. In 
such a case it will likely be expedient to simply set /3 = 0. Further discussion of these issues 
may be found below. 

With the choice of trial density fflTD . the Lagrangian is easily evaluated up to an irrelevant 
additive constant as 



(16) 


9ij — {QijQj) 

Mi = {L*Qi) 

(^ = A' {L*QiL*Q,) 
X, = {{Q,-{Q,))T) 
Yi = {{L*Qi)T) 



The non-negative dehnite matrix g{X) is the Fisher information matrix for the random vector 
Q whose trial densities are the given maximum entropy statistical model family. Thus it is 
also the metric tensor of the trial density manifold. Note also that the scalar function F here 
will vanish identically if ip happens to be chosen an invariant function for the dynamical 
system. In the case it is a quadratic approximation to such an invariant then F will be 
expected to be small and hence the linear terms in the Lagrangian will be also. Finally 
it is worth emphasizing that due to flT^ . the action is always bounded below and so the 
Euclidean path integral we have proposed here is well dehned despite the complexity of 
equation ([15]). 
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It is clear now that our first and second requirements for a path measure from section 2 
above have been met since all these fields are specified as functions of A (and hence the mo¬ 
ments of the slow variables via a Legendre transformation) and a Lagrangian defined based 
on a clear information theoretic formulation which generalizes entropy from the equilibrium 
Boltzmann principle. The third requirement holds in a certain sense when the system is 
Hamiltonian but the more general autonomous case considered here requires further inves¬ 
tigation. 


B. Paths and endpoints 


We have seen how to associate a weight to a path using the information loss action. It 
is clear however that we further require a non-negative weight for any trial density at a 
prescribed endpoint time in order to evaluate it’s significance as an approximation for the 
true probability density at that time. The path integral proposed by Onsager and Machlup 
and generalized here, clearly achieves this objective in a natural manner (see equation (1^ ). 
The trial density weight at a fixed time is termed a consistency distribution and given it’s 
formulation as a path integral is analogous to a quantum wavefunction. Indeed it satisfies a 
(Wick rotated) Schrodinger equation (see jl|). 

It is important to emphasize the significance of the path integral here. It is the sum of all 
path weights with fixed endpoints that is used to define the consistency distribution. Such a 
path integral is not in general a monotonic function of the extremal action which applies 
to a path satisfying the Euler-Lagrange equations for Such an action is termed the 
classical action Sd in the path integral literature. In the special case that the Lagrangian 
.if is a quadratic function of A it may however be proven (see e.g. {3) Chapter 6) that 


J DXW [X] = Bexp{-Sd) 


with B not dependent on the value of A at the path endpoint. Clearly in that simple case 
(which also happens to be the original Onsager-Machlup one) only the extremal actions are 
important in defining the consistency distribution. In most realistic situations however .if 
is not quadratic and so this simplifying situation may not apply. 

More concretely, one can imagine a scenario when the extremal paths between two sets 
of endpoints have the same minimal action but the path integral is quite different. This 
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situation is sketched in Figure 1. 



Figure 1. The importance of a path integral versus extremal actions in defining the consistency 
distribution. In the example sketch we have the situation that IF 3 , W 4 << Wi, W 2 < Wmax which 
is when “quantum” effects are significant. 


V. DISCUSSION 

The present contribution is intended as a documentation of a closure formalism which 
works in principle for a wide class of realistic dynamical systems including those pertain¬ 
ing to turbulence. We discuss now the practical problems involved in implementing this 
formalism. 

One issue to consider is the choice of approximating or trial densities. In section 4 we 
posited that some kind of exponential family of densities was appropriate and this assump¬ 
tion allowed the formal completion of the calculation of the Lagrangian. A key consideration 
here is whether a member of the chosen trial density manifold could be considered to be a 
reasonable approximation to the non-equilibrium densities of the full system. If this is not 
so one should not expect a priori the present closure to perform well. This is an important 
issue since a complete calculation of the Lagrangian involves knowledge of the trial density 
partition function Z as well as its various moments. In the case that a Gaussian trial density 
is a good approximating density then this calculation is straightforward. It is clear however 
that a more general class of trial densities may need to be chosen carefully in order to bal- 
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ance approximating accuracy against the ability to analytically evaluate moments and the 
partition function. 

It is also important to consider the statistical behaviour of the unresolved or fast variables 
within the system. The implicit assumption of the present closure method is that on the slow 
timescale of interest, the fast variable statistics of the system are effectively relaxed to those 
applicable for an equilibrium density. Referring to the ansatz trial density (fTT|) the only 
dependency of the density on fast variables is assumed to come through the function ■^(x). 
For a forced dissipative turbulent system the obvious choice for this function is one such that 
C'exp(—/3'0(x)) is close to the equilibrium density of the system since that ensures the trial 
densities have the correct statistical properties for the fast variables (at least on the timescale 
of interest). Knowledge of such a density may be in general difficult to obtain theoretically 
and instead require a direct numerical simulation of the system. For a decaying turbulent 
system, the asymptotic state is usually one of no motion in which case the exp{—/3'ip{x)) 
piece could plausibly be replaced simply by the products of delta functions for the fast 
variables since these are the limits of Gaussian densities as the variance go to zero. 


A further issue to consider is the numerical evaluation of the path integrals proposed 
here. General Markov Ghain Monte Garlo methods for such integrals are well known (see, 
and these have been very recently adapted by the author to the present 
15|). One constraint here is that the number of coarse graining variables used 


for example, Hi 
situation (see 


needs to be manageably small in order that the numerical calculations remain tractable. A 
further constraint concerns the time discretisation used to approximate the action integral 
since the numerical cost of the MGMG method depends significantly on the number of time 
ordinates used as documented by Geperley and co-workers. It should be noted that a poten¬ 
tial advantage of the path integral approach as opposed to other methods used previously by 
the author and collaborator (see j^) is that it is applicable to far from statistical equilibrium 
situations. The other method relies on a Taylor expansion of an associated Hamilton-Jacobi 
equation about the equilibrium density. 
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